
EnableExplicit

;{ other
	Macro Math_DoubleQuotes
		"
	EndMacro
;}
;{ Structures
	Structure vec2
		x.f
		y.f
	EndStructure
	Structure vec3
		x.f
		y.f
		z.f
	EndStructure
	Structure vec4
		x.f
		y.f
		z.f
		w.f
	EndStructure
	Structure mat2_
		_11.f : _12.f
		_21.f : _22.f
	EndStructure
	Structure mat3
		_11.f : _12.f : _13.f
		_21.f : _22.f : _23.f
		_31.f : _32.f : _33.f
	EndStructure
	Structure mat4
		_11.f : _12.f : _13.f : _14.f
		_21.f : _22.f : _23.f : _24.f
		_31.f : _32.f : _33.f : _34.f
		_41.f : _42.f : _43.f : _44.f
	EndStructure
;}

;{ Integers
	Procedure.l sInt_Random(one.l, two.l)
		If one <> two
			If one < two
				ProcedureReturn Random(two-one) + one
			Else
				ProcedureReturn Random(one-two) + two
			EndIf
		Else
			ProcedureReturn one
		EndIf
	EndProcedure
;}
;{ Reals
	Procedure.f Real_Random(one.f, two.f)
		ProcedureReturn sInt_Random(one*1000000, two*1000000)*0.000001
	EndProcedure
;}
;{ Vectors
	;{ variables
		Threaded _vec_normalize_k.f
		Global Vec2_UnityX.vec4 : Vec2_UnityX\x = 1
		Global Vec2_UnityY.vec4 : Vec2_UnityY\y = 1
		Global Vec3_UnityX.vec4 : Vec3_UnityX\x = 1
		Global Vec3_UnityY.vec4 : Vec3_UnityY\y = 1
		Global Vec3_UnityZ.vec4 : Vec3_UnityZ\z = 1
		Global Vec4_UnityX.vec4 : Vec4_UnityX\x = 1
		Global Vec4_UnityY.vec4 : Vec4_UnityY\y = 1
		Global Vec4_UnityZ.vec4 : Vec4_UnityZ\z = 1
		Global Vec4_UnityW.vec4 : Vec4_UnityW\w = 1
	;}
	;{ 2d
		Macro				Vec2_Set(a, ax, ay)
			a\x = ax
			a\y = ay
		EndMacro
		Macro				Vec2_Copy(t, s)
			t\x = s\x
			t\y = s\y
		EndMacro
		Macro				Vec2_Scale(t, s, k)
			t\x = s\x * k
			t\y = s\y * k
		EndMacro
		Macro				Vec2_Add(t, s1, s2)
			t\x = s1\x + s2\x
			t\y = s1\y + s2\y
		EndMacro
		Macro				Vec2_Sub(t, s1, s2)
			t\x = s1\x - s2\x
			t\y = s1\y - s2\y
		EndMacro
		Macro				Vec2_Mul(c, a, b)
			c\x = a\x * b\x
			c\y = a\y * b\y
		EndMacro
		Macro				Vec2_Dot(a, b)
			(a\x*b\x + a\y*b\y)
		EndMacro
		Macro				Vec2_LenQ(a)
			Vec2_Dot(a, a)
		EndMacro
		Macro				Vec2_Len(a)
			Sqr(Vec2_LenQ(a))
		EndMacro
		Macro				Vec2_Normalize(a)
			_vec_normalize_k = 1/Vec2_Len(a)
			Vec2_Scale(a, a, _vec_normalize_k)
		EndMacro
		Macro				Vec2_Orth(t, s, k)
			t\x = -s\y*k
			t\y = s\x*k
		EndMacro
		Macro				Vec2_Debug(v)
			Debug Math_DoubleQuotes#v#Math_DoubleQuotes + ": [" + StrF(v\x, 6) + ", " + StrF(v\y, 6) + "]"
		EndMacro
	;}
	;{ 3d
		Macro				Vec3_Set(a, ax, ay, az)
			a\x = ax
			a\y = ay
			a\z = az
		EndMacro
		Macro				Vec3_Random(t)
			t\x = Real_Random(-1, 1)
			t\y = Real_Random(-1, 1)
			t\z = Real_Random(-1, 1)
		EndMacro
		Macro				Vec3_Vec(t, s, _z)
			t\x = s\x
			t\y = s\y
			t\z = _z
		EndMacro
		Macro				Vec3_Copy(t, s)
			t\x = s\x
			t\y = s\y
			t\z = s\z
		EndMacro
		Macro				Vec3_Scale(t, s, k)
			t\x = s\x * k
			t\y = s\y * k
			t\z = s\z * k
		EndMacro
		Macro				Vec3_Add(t, s1, s2)
			t\x = s1\x + s2\x
			t\y = s1\y + s2\y
			t\z = s1\z + s2\z
		EndMacro
		Macro				Vec3_Sub(t, s1, s2)
			t\x = s1\x - s2\x
			t\y = s1\y - s2\y
			t\z = s1\z - s2\z
		EndMacro
		Macro				Vec3_Mul(c, a, b)
			c\x = a\x * b\x
			c\y = a\y * b\y
			c\z = a\z * b\z
		EndMacro
		Macro				Vec3_Dot(a, b)
			((a\x * b\x) + (a\y * b\y) + (a\z * b\z))
		EndMacro
		Macro				Vec3_LenQ(a)
			Vec3_Dot(a, a)
		EndMacro
		Macro				Vec3_Len(a)
			Sqr(Vec3_LenQ(a))
		EndMacro
		Macro				Vec3_Normalize(a)
			_vec_normalize_k = 1/Vec3_Len(a)
			Vec3_Scale(a, a, _vec_normalize_k)
		EndMacro
		Macro				Vec3_Cross(c, a, b)
			c\x = a\y*b\z - a\z*b\y
			c\y =  a\z*b\x - a\x*b\z
			c\z =  a\x*b\y - a\y*b\x
		EndMacro
		Macro				Vec3_Debug(v)
			Debug Math_DoubleQuotes#v#Math_DoubleQuotes + ": [" + StrF(v\x) + ", " + StrF(v\y) + ", " +StrF(v\z) + "]"
		EndMacro
	;}
	;{ 4d
		Macro				Vec4_Set(a, ax, ay, az, aw)
			a\x = ax
			a\y = ay
			a\z = az
			a\w  = aw
		EndMacro
		Macro				Vec4_Vec(t, s, _w)
			t\x = s\x
			t\y = s\y
			t\z = s\z
			t\w = _w
		EndMacro
		Macro				Vec4_Copy(t, s)
			t\x = s\x
			t\y = s\y
			t\z = s\z
			t\w = s\w
		EndMacro
		Macro				Vec4_Scale(t, s, k)
			t\x = s\x * k
			t\y = s\y * k
			t\z = s\z * k
			t\w = s\w * k
		EndMacro
		Macro				Vec4_Add(t, s1, s2)
			t\x = s1\x + s2\x
			t\y = s1\y + s2\y
			t\z = s1\z + s2\z
			t\w = s1\w + s2\w
		EndMacro
		Macro				Vec4_Sub(t, s1, s2)
			t\x = s1\x - s2\x
			t\y = s1\y - s2\y
			t\z = s1\z - s2\z
			t\w = s1\w - s2\w
		EndMacro
		Macro				Vec4_Mul(c, a, b)
			c\x = a\x * b\x
			c\y = a\y * b\y
			c\z = a\z * b\z
			c\w = a\w * b\w
		EndMacro
		Macro				Vec4_Dot(a, b)
			((a\x*b\x) + (a\y*b\y) + (a\z*b\z) + (a\w*b\w))
		EndMacro
		Macro				Vec4_LenQ(a)
			Vec4_Dot(a, a)
		EndMacro
		Macro				Vec4_Len(a)
			Sqr(Vec4_LenQ(a))
		EndMacro
		Macro				Vec4_Normalize(a)
			_vec_normalize_k = 1/Vec4_Len(a)
			Vec4_Scale(a, a, _vec_normalize_k)
		EndMacro
		Macro				Vec4_Hamilton(c, a, b)
			c\x = a\w * b\x + a\x * b\w + a\y * b\z - a\z * b\y
			c\y = a\w * b\y - a\x * b\z + a\y * b\w + a\z * b\x
			c\z = a\w * b\z + a\x * b\y - a\y * b\x + a\z * b\w
			c\w = a\w * b\w - a\x * b\x - a\y * b\y - a\z * b\z
		EndMacro
		Macro				Vec4_Debug(v)
			Debug Math_DoubleQuotes#v#Math_DoubleQuotes + ": [" + StrF(v\x, 6) + ", " + StrF(v\y, 6) + ", " +StrF(v\z, 6) + ", " + StrF(v\w, 6) + "]"
		EndMacro
	;}
;}
;{ Matrices
	;{ variables
		Threaded _mat3_mulm_Z.mat3
		Threaded _mat3_mulv_c.vec3
		Threaded _mat4_mulm_Z.mat4
	;}
	;{ 2x2
		Macro			Mat2_Set(A, _11_, _12_, _21_, _22_)
			A\_11 = _11_
			A\_12 = _12_
			A\_21 = _21_
			A\_22 = _22_
		EndMacro
		Macro			Mat2_MulM(Z, X, Y)
			Z\_11 = (X\_11 * Y\_11)  +  (X\_12 * Y\_21)
			Z\_12 = (X\_11 * Y\_12)  +  (X\_12 * Y\_22)
			Z\_21 = (X\_21 * Y\_11)  +  (X\_22 * Y\_21)
			Z\_22 = (X\_21 * Y\_12)  +  (X\_22 * Y\_22)
		EndMacro
		Macro			Mat2_MulV(c, A, b)
			c\x = A\_11 * b\x  +  A\_12 * b\y
			c\y = A\_21 * b\x  +  A\_22 * b\y
		EndMacro
		Macro			Mat2_Scale(T, S, k)
			T\_11 = S\_11 * k
			T\_12 = S\_12 * k
			T\_21 = S\_21 * k
			T\_22 = S\_22 * k
		EndMacro
		Macro			Mat2_Transpose(T, S)
			T\_11 = S\_11
			T\_12 = S\_21
			T\_21 = S\_12
			T\_22 = S\_22
		EndMacro
		Macro			Mat2_CreateIdentity(A)
			Mat2_Set(A, 1, 0, 0, 1)
		EndMacro
		Procedure	Mat2_CreateRotation(*M.mat2_, a.f)
			Protected.f cos, sin
			cos = Cos(a)
			sin = Sin(a)
			*M\_11 = cos
			*M\_12 = -sin
			*M\_21 = sin
			*M\_22 = cos
		EndProcedure
		Procedure	Mat2_CreateScaleComp(*M.mat2_, kx.f, ky.f)
			*M\_11 = kx
			*M\_12 = 0
			*M\_21 = 0
			*M\_22 = ky
		EndProcedure
		Procedure	Mat2_CreateScale(*M.mat2_, *scl.vec2)
			Mat2_CreateScaleComp(*M, *scl\x, *scl\y)
		EndProcedure
		Macro 			Mat2_Debug(A)
			Debug StrF(A\_11, 1) + ", " + StrF(A\_12, 1)
			Debug StrF(A\_21, 1) + ", " + StrF(A\_22, 1)
		EndMacro
	;}
	;{ 3x3
		Macro			Mat3_SetL1(A, _1, _2, _3)
			A\_11 = _1
			A\_12 = _2
			A\_13 = _3
		EndMacro
		Macro			Mat3_SetL2(A, _1, _2, _3)
			A\_21 = _1
			A\_22 = _2
			A\_23 = _3
		EndMacro
		Macro			Mat3_SetL3(A, _1, _2, _3)
			A\_31 = _1
			A\_32 = _2
			A\_33 = _3
		EndMacro
		Procedure	Mat3_SetV(*M.mat3, *a.vec3, *b.vec3, *c.vec3)
			Mat3_SetL1(*M,  *a\x,  *b\x, *c\x)
			Mat3_SetL2(*M,  *a\y,  *b\y, *c\y)
			Mat3_SetL3(*M,  *a\z,  *b\z,  *c\z)
		EndProcedure
		Macro			Mat3_Copy(T, S)
			T\_11 = S\_11
			T\_12 = S\_12
			T\_13 = S\_13
			T\_21 = S\_21
			T\_22 = S\_22
			T\_23 = S\_23
			T\_31 = S\_31
			T\_32 = S\_32
			T\_33 = S\_33
		EndMacro
		Procedure	Mat3_MulM(*Z.mat3, *X.mat3, *Y.mat3)
			If Not(*Z = *X And *Z = *Y)
				*Z\_11 = *X\_11 * *Y\_11 + *X\_12 * *Y\_21 + *X\_13 * *Y\_31
				*Z\_12 = *X\_11 * *Y\_12 + *X\_12 * *Y\_22 + *X\_13 * *Y\_32
				*Z\_13 = *X\_11 * *Y\_13 + *X\_12 * *Y\_23 + *X\_13 * *Y\_33
				*Z\_21 = *X\_21 * *Y\_11 + *X\_22 * *Y\_21 + *X\_23 * *Y\_31
				*Z\_22 = *X\_21 * *Y\_12 + *X\_22 * *Y\_22 + *X\_23 * *Y\_32
				*Z\_23 = *X\_21 * *Y\_13 + *X\_22 * *Y\_23 + *X\_23 * *Y\_33
				*Z\_31 = *X\_31 * *Y\_11 + *X\_32 * *Y\_21 + *X\_33 * *Y\_31
				*Z\_32 = *X\_31 * *Y\_12 + *X\_32 * *Y\_22 + *X\_33 * *Y\_32
				*Z\_33 = *X\_31 * *Y\_13 + *X\_32 * *Y\_23 + *X\_33 * *Y\_33
			Else
				If Not(*Z = _mat3_mulm_Z)
					Mat3_MulM(_mat3_mulm_Z, *X, *Y)
					Mat3_Copy(*Z, _mat3_mulm_Z)
				Else
					Debug Str(#PB_Compiler_Line) + " line: recursive call here"
				EndIf
			EndIf
		EndProcedure
		Procedure	Mat3_MulV(*c.vec3, *A.mat3, *b.vec3)
			If Not *c = *b
				*c\x = *A\_11 * *b\x  +  *A\_12 * *b\y  +  *A\_13 * *b\z
				*c\y = *A\_21 * *b\x  +  *A\_22 * *b\y  +  *A\_23 * *b\z
				*c\z = *A\_31 * *b\x  +  *A\_32 * *b\y  +  *A\_33 * *b\z
			Else
				If Not(*c = _mat3_mulv_c)
					Mat3_MulV(_mat3_mulv_c, *A, *b)
					Vec3_Copy(*c, _mat3_mulv_c)
				Else
					Debug Str(#PB_Compiler_Line) + " line: recursive call here"
				EndIf
			EndIf
		EndProcedure
		Procedure	Mat3_Scale(*T.mat3, *S.mat3, k.f)
			*T\_11 = *S\_11 * k
			*T\_12 = *S\_12 * k
			*T\_13 = *S\_13 * k
			*T\_21 = *S\_21 * k
			*T\_22 = *S\_22 * k
			*T\_23 = *S\_23 * k
			*T\_31 = *S\_31 * k
			*T\_32 = *S\_32 * k
			*T\_33 = *S\_33 * k
		EndProcedure
		Procedure	Mat3_Transpose(*T.mat3, *S.mat3)
			*T\_11 = *S\_11
			*T\_12 = *S\_21
			*T\_13 = *S\_31
			*T\_21 = *S\_12
			*T\_22 = *S\_22
			*T\_23 = *S\_32
			*T\_31 = *S\_13
			*T\_32 = *S\_23
			*T\_33 = *S\_33
		EndProcedure
		Macro			Mat3_Det(M)
			(M\_11*(M\_22*M\_33 - M\_32*M\_23) - M\_12*(M\_21*M\_33 - M\_31*M\_23) + M\_13*(M\_21*M\_32 - M\_31*M\_22))
		EndMacro
		Procedure	Mat3_Inverse(*T.mat3, *S.mat3)
			Protected.f k_det
			k_det = 1/Mat3_Det(*S)
			*T\_11 =    (*S\_22 * *S\_33 - *S\_23 * *S\_32) * k_det
			*T\_12 =  -(*S\_12 * *S\_33 - *S\_32 * *S\_13) * k_det
			*T\_13 =    (*S\_12 * *S\_23 - *S\_22 * *S\_13) * k_det
			
			*T\_21 =  -(*S\_21 * *S\_33 - *S\_23 * *S\_31) * k_det
			*T\_22 =    (*S\_11 * *S\_33 - *S\_31 * *S\_13) * k_det
			*T\_23 =  -(*S\_11 * *S\_23 - *S\_21 * *S\_13) * k_det
			
			*T\_31 =    (*S\_21 * *S\_32 - *S\_31 * *S\_22) * k_det
			*T\_32 =  -(*S\_11 * *S\_32 - *S\_31 * *S\_12) * k_det
			*T\_33 =    (*S\_11 * *S\_22 - *S\_12 * *S\_21) * k_det
		EndProcedure
		Macro			Mat3_CreateIdentity(M)
			Mat3_SetL1(M, 1, 0, 0)
			Mat3_SetL2(M, 0, 1, 0)
			Mat3_SetL3(M, 0, 0, 1)
		EndMacro
		Procedure	Mat3_CreateRandom(*M.mat3)
			Mat3_SetL1(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
			Mat3_SetL2(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
			Mat3_SetL3(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
		EndProcedure
		Procedure	Mat3_CreateRotationZ(*M.mat3, a.f)
			Mat3_SetL1(*M, Cos(a), -Sin(a), 0)
			Mat3_SetL2(*M, Sin(a),   Cos(a), 0)
			Mat3_SetL3(*M,     0     ,       0      , 1)
		EndProcedure
		Procedure	Mat3_CreateRotationY(*M.mat3, a.f)
			; fixe - sign may be wrong
			Mat3_SetL1(*M, Cos(a), 0, Sin(a))
			Mat3_SetL2(*M,       0     , 1,      0    )
			Mat3_SetL3(*M,-Sin(a),  0, Cos(a))
		EndProcedure
		Procedure	Mat3_CreateRotationX(*M.mat3, a.f)
			Mat3_SetL1(*M,	1,      0     ,       0     )
			Mat3_SetL2(*M,	0, Cos(a), -Sin(a))
			Mat3_SetL3(*M,	0, Sin(a), Cos(a))
		EndProcedure
		Procedure	Mat3_CreateRotationComp(*M.mat3, rx.f, ry.f, rz.f)
			Protected.mat3 X, Y, Z, YX
			Mat3_CreateRotationX(X, rx)
			Mat3_CreateRotationY(Y, ry)
			Mat3_CreateRotationZ(Z, rz)
			Mat3_MulM(YX, Y, X)
			Mat3_MulM(*M, Z, YX)
		EndProcedure
		Procedure	Mat3_CreateRotation(*M.mat3, *r.vec3)
			Mat3_CreateRotationComp(*M, *r\x, *r\y, *r\z)
		EndProcedure
		Procedure	Mat3_CreateScaleComp(*M.mat3, kx.f, ky.f, kz.f)
			Mat3_SetL1(*M, kx, 0, 0)
			Mat3_SetL2(*M, 0, ky, 0)
			Mat3_SetL3(*M, 0, 0, kz)
		EndProcedure
		Procedure	Mat3_CreateScale(*M, *k.vec3)
			Mat3_CreateScaleComp(*M, *k\x, *k\y, *k\z)
		EndProcedure
		Procedure	Mat3_FromQuat(*M.mat3, *Q.vec4)
			Protected.f x, y, z, w, zz2, yy2, xx2, xy2, zw2, xz2, yw2, yz2, xw2
			x = *Q\x : y = *Q\y : z = *Q\z : w = *Q\w
			xx2 = x*x*2
			xy2 = x*y*2
			xz2 = x*z*2
			xw2 = x*w*2
			yy2 = y*y*2
			yz2 = y*z*2
			yw2 = y*w*2
			zz2 = z*z*2
			zw2 = z*w*2
			Mat3_SetL1(*M, 1 - yy2 - zz2, xy2 - zw2, xz2 + yw2)
			Mat3_SetL2(*M, xy2 + zw2, 1 - xx2 - zz2, yz2 - xw2)
			Mat3_SetL3(*M, xz2 - yw2, yz2 + xw2, 1 - xx2 - yy2)
		EndProcedure
		Procedure	Mat3_Debug(*A.mat3)
			Debug StrF(*A\_11, 1) + ", " + StrF(*A\_12, 1) + ", " + StrF(*A\_13, 1)
			Debug StrF(*A\_21, 1) + ", " + StrF(*A\_22, 1) + ", " + StrF(*A\_23, 1)
			Debug StrF(*A\_31, 1) + ", " + StrF(*A\_32, 1) + ", " + StrF(*A\_33, 1)
		EndProcedure
	;}
	;{ 4x4
		Macro 				Mat4_SetL1(A, _1, _2, _3, _4)
			A\_11 = _1
			A\_12 = _2
			A\_13 = _3
			A\_14 = _4
		EndMacro
		Macro 				Mat4_SetL2(A, _1, _2, _3, _4)
			A\_21 = _1
			A\_22 = _2
			A\_23 = _3
			A\_24 = _4
		EndMacro
		Macro 				Mat4_SetL3(A, _1, _2, _3, _4)
			A\_31 = _1
			A\_32 = _2
			A\_33 = _3
			A\_34 = _4
		EndMacro
		Macro 				Mat4_SetL4(A, _1, _2, _3, _4)
			A\_41 = _1
			A\_42 = _2
			A\_43 = _3
			A\_44 = _4
		EndMacro
		Procedure		Mat4_SetV(*M.mat4, *A.vec4, *B.vec4, *C.vec4, *D.vec4)
			Mat4_SetL1(*M,  *A\x,  *B\x, *C\x,  *D\x)
			Mat4_SetL2(*M,  *A\y,  *B\y, *C\y,  *D\y)
			Mat4_SetL3(*M,  *A\z,  *B\z,  *C\z,  *D\z)
			Mat4_SetL4(*M, *A\w, *B\w, *C\w, *D\w)
		EndProcedure
		Procedure		Mat4_SetM(*C.mat4, *A.mat3, *b.vec4)
			Mat4_SetL1(*C, *A\_11, *A\_12, *A\_13, *b\x)
			Mat4_SetL2(*C, *A\_21, *A\_22, *A\_23, *b\y)
			Mat4_SetL3(*C, *A\_31, *A\_32, *A\_33, *b\z)
			Mat4_SetL4(*C,       0	  ,	   0	  ,	   0	  , *b\w)
		EndProcedure
		Macro				Mat4_Copy(T, S)
			T\_11 = S\_11
			T\_12 = S\_12
			T\_13 = S\_13
			T\_14 = S\_14
			T\_21 = S\_21
			T\_22 = S\_22
			T\_23 = S\_23
			T\_24 = S\_24
			T\_31 = S\_31
			T\_32 = S\_32
			T\_33 = S\_33
			T\_34 = S\_34
			T\_41 = S\_41
			T\_42 = S\_42
			T\_43 = S\_43
			T\_44 = S\_44
		EndMacro
		Procedure		Mat4_MulM(*Z.mat4, *X.mat4, *Y.mat4)
			If Not(*Z = *X And *Z = *Y)
				*Z\_11 = *X\_11 * *Y\_11  +  *X\_12 * *Y\_21  +  *X\_13 * *Y\_31  +  *X\_14 * *Y\_41
				*Z\_12 = *X\_11 * *Y\_12  +  *X\_12 * *Y\_22  +  *X\_13 * *Y\_32  +  *X\_14 * *Y\_42
				*Z\_13 = *X\_11 * *Y\_13  +  *X\_12 * *Y\_23  +  *X\_13 * *Y\_33  +  *X\_14 * *Y\_43
				*Z\_14 = *X\_11 * *Y\_14  +  *X\_12 * *Y\_24  +  *X\_13 * *Y\_34  +  *X\_14 * *Y\_44
				*Z\_21 = *X\_21 * *Y\_11  +  *X\_22 * *Y\_21  +  *X\_23 * *Y\_31  +  *X\_24 * *Y\_41
				*Z\_22 = *X\_21 * *Y\_12  +  *X\_22 * *Y\_22  +  *X\_23 * *Y\_32  +  *X\_24 * *Y\_42
				*Z\_23 = *X\_21 * *Y\_13  +  *X\_22 * *Y\_23  +  *X\_23 * *Y\_33  +  *X\_24 * *Y\_43
				*Z\_24 = *X\_21 * *Y\_14  +  *X\_22 * *Y\_24  +  *X\_23 * *Y\_34  +  *X\_24 * *Y\_44
				*Z\_31 = *X\_31 * *Y\_11  +  *X\_32 * *Y\_21  +  *X\_33 * *Y\_31  +  *X\_34 * *Y\_41
				*Z\_32 = *X\_31 * *Y\_12  +  *X\_32 * *Y\_22  +  *X\_33 * *Y\_32  +  *X\_34 * *Y\_42
				*Z\_33 = *X\_31 * *Y\_13  +  *X\_32 * *Y\_23  +  *X\_33 * *Y\_33  +  *X\_34 * *Y\_43
				*Z\_34 = *X\_31 * *Y\_14  +  *X\_32 * *Y\_24  +  *X\_33 * *Y\_34  +  *X\_34 * *Y\_44
				*Z\_41 = *X\_41 * *Y\_11  +  *X\_42 * *Y\_21  +  *X\_43 * *Y\_31  +  *X\_44 * *Y\_41
				*Z\_42 = *X\_41 * *Y\_12  +  *X\_42 * *Y\_22  +  *X\_43 * *Y\_32  +  *X\_44 * *Y\_42
				*Z\_43 = *X\_41 * *Y\_13  +  *X\_42 * *Y\_23  +  *X\_43 * *Y\_33  +  *X\_44 * *Y\_43
				*Z\_44 = *X\_41 * *Y\_14  +  *X\_42 * *Y\_24  +  *X\_43 * *Y\_34  +  *X\_44 * *Y\_44
			Else
				If Not(*Z = _mat4_mulm_Z)
					Mat4_MulM(_mat4_mulm_Z, *X, *Y)
					Mat4_Copy(*Z, _mat4_mulm_Z)
				Else
					Debug Str(#PB_Compiler_Line) + " line: recursive call here"
				EndIf
			EndIf
		EndProcedure
		Procedure		Mat4_MulV(*C.vec4, *A.mat4, *B.vec4)
			*C\x = *A\_11 * *B\x  +  *A\_12 * *B\y  +  *A\_13 * *B\z  +  *A\_14 * *B\w
			*C\y = *A\_21 * *B\x  +  *A\_22 * *B\y  +  *A\_23 * *B\z  +  *A\_24 * *B\w
			*C\z = *A\_31 * *B\x  +  *A\_32 * *B\y  +  *A\_33 * *B\z  +  *A\_34 * *B\w
			*C\z = *A\_41 * *B\x  +  *A\_42 * *B\y  +  *A\_43 * *B\z  +  *A\_44 * *B\w
		EndProcedure
		Procedure		Mat4_Scale(*T.mat4, *S.mat4, k.f)
			*T\_11 = *S\_11 * k
			*T\_12 = *S\_12 * k
			*T\_13 = *S\_13 * k
			*T\_14 = *S\_14 * k
			*T\_21 = *S\_21 * k
			*T\_22 = *S\_22 * k
			*T\_23 = *S\_23 * k
			*T\_24 = *S\_24 * k
			*T\_31 = *S\_31 * k
			*T\_32 = *S\_32 * k
			*T\_33 = *S\_33 * k
			*T\_34 = *S\_34 * k
			*T\_41 = *S\_41 * k
			*T\_42 = *S\_42 * k
			*T\_43 = *S\_43 * k
			*T\_44 = *S\_44 * k
		EndProcedure
		Procedure		Mat4_Transpose(*T.mat4, *S.mat4)
			*T\_11 = *S\_11
			*T\_12 = *S\_21
			*T\_13 = *S\_31
			*T\_14 = *S\_41
			*T\_21 = *S\_12
			*T\_22 = *S\_22
			*T\_23 = *S\_32
			*T\_24 = *S\_42
			*T\_31 = *S\_13
			*T\_32 = *S\_23
			*T\_33 = *S\_33
			*T\_34 = *S\_43
			*T\_41 = *S\_14
			*T\_42 = *S\_24
			*T\_43 = *S\_34
			*T\_44 = *S\_44
		EndProcedure
		Procedure.f	Mat4_Minor(*M.mat4, i.l, j.l)
			Protected.mat3 m
			If i = 1
				If j = 1
					Mat3_SetL1(m, *M\_22, *M\_23, *M\_24)
					Mat3_SetL2(m, *M\_32, *M\_33, *M\_34)
					Mat3_SetL3(m, *M\_42, *M\_43, *M\_44)
				EndIf
				If j = 2
					Mat3_SetL1(m, *M\_21, *M\_23, *M\_24)
					Mat3_SetL2(m, *M\_31, *M\_33, *M\_34)
					Mat3_SetL3(m, *M\_41, *M\_43, *M\_44)
				EndIf
				If j = 3
					Mat3_SetL1(m, *M\_21, *M\_22, *M\_24)
					Mat3_SetL2(m, *M\_31, *M\_32, *M\_34)
					Mat3_SetL3(m, *M\_41, *M\_42, *M\_44)
				EndIf
				If j = 4
					Mat3_SetL1(m, *M\_21, *M\_22, *M\_23)
					Mat3_SetL2(m, *M\_31, *M\_32, *M\_33)
					Mat3_SetL3(m, *M\_41, *M\_42, *M\_43)
				EndIf
			EndIf
			If i = 2
				If j = 1
					Mat3_SetL1(m, *M\_12, *M\_13, *M\_14)
					Mat3_SetL2(m, *M\_32, *M\_33, *M\_34)
					Mat3_SetL3(m, *M\_42, *M\_43, *M\_44)
				EndIf
				If j = 2
					Mat3_SetL1(m, *M\_11, *M\_13, *M\_14)
					Mat3_SetL2(m, *M\_31, *M\_33, *M\_34)
					Mat3_SetL3(m, *M\_41, *M\_43, *M\_44)
				EndIf
				If j = 3
					Mat3_SetL1(m, *M\_11, *M\_12, *M\_14)
					Mat3_SetL2(m, *M\_31, *M\_32, *M\_34)
					Mat3_SetL3(m, *M\_41, *M\_42, *M\_44)
				EndIf
				If j = 4
					Mat3_SetL1(m, *M\_11, *M\_12, *M\_13)
					Mat3_SetL2(m, *M\_31, *M\_32, *M\_33)
					Mat3_SetL3(m, *M\_41, *M\_42, *M\_43)
				EndIf
			EndIf
			If i = 3
				If j = 1
					Mat3_SetL1(m, *M\_12, *M\_13, *M\_14)
					Mat3_SetL2(m, *M\_22, *M\_23, *M\_24)
					Mat3_SetL3(m, *M\_42, *M\_43, *M\_44)
				EndIf
				If j = 2
					Mat3_SetL1(m, *M\_11, *M\_13, *M\_14)
					Mat3_SetL2(m, *M\_21, *M\_23, *M\_24)
					Mat3_SetL3(m, *M\_41, *M\_43, *M\_44)
				EndIf
				If j = 3
					Mat3_SetL1(m, *M\_11, *M\_12, *M\_14)
					Mat3_SetL2(m, *M\_21, *M\_22, *M\_24)
					Mat3_SetL3(m, *M\_41, *M\_42, *M\_44)
				EndIf
				If j = 4
					Mat3_SetL1(m, *M\_11, *M\_12, *M\_13)
					Mat3_SetL2(m, *M\_21, *M\_22, *M\_23)
					Mat3_SetL3(m, *M\_41, *M\_42, *M\_43)
				EndIf
			EndIf
			If i = 4
				If j = 1
					Mat3_SetL1(m, *M\_12, *M\_13, *M\_14)
					Mat3_SetL2(m, *M\_22, *M\_23, *M\_24)
					Mat3_SetL3(m, *M\_32, *M\_33, *M\_34)
				EndIf
				If j = 2
					Mat3_SetL1(m, *M\_11, *M\_13, *M\_14)
					Mat3_SetL2(m, *M\_21, *M\_23, *M\_24)
					Mat3_SetL3(m, *M\_31, *M\_33, *M\_34)
				EndIf
				If j = 3
					Mat3_SetL1(m, *M\_11, *M\_12, *M\_14)
					Mat3_SetL2(m, *M\_21, *M\_22, *M\_24)
					Mat3_SetL3(m, *M\_31, *M\_32, *M\_34)
				EndIf
				If j = 4
					Mat3_SetL1(m, *M\_11, *M\_12, *M\_13)
					Mat3_SetL2(m, *M\_21, *M\_22, *M\_23)
					Mat3_SetL3(m, *M\_31, *M\_32, *M\_33)
				EndIf
			EndIf
			ProcedureReturn Mat3_Det(m)
		EndProcedure
		Procedure.f	Mat4_Det(*M.mat4)
			Protected.f a1, a2, a3, a4
			a1 =    *M\_11*Mat4_Minor(*M, 1, 1)
			a2 = -(*M\_12*Mat4_Minor(*M, 1, 2))
			a3 =    *M\_13*Mat4_Minor(*M, 1, 3)
			a4 = -(*M\_14*Mat4_Minor(*M, 1, 4))
			ProcedureReturn a1 + a2 + a3 + a4
		EndProcedure
		Procedure		Mat4_Inverse(*T.mat4, *S.mat4)
			Protected.mat4 A
			Protected.f det, k_det
			A\_11 =   Mat4_Minor(*S, 1, 1)
			A\_12 =  -Mat4_Minor(*S, 1, 2)
			A\_13 =   Mat4_Minor(*S, 1, 3)
			A\_14 = -Mat4_Minor(*S, 1, 4)
			A\_21 = -Mat4_Minor(*S, 2, 1)
			A\_22 =  Mat4_Minor(*S, 2, 2)
			A\_23 = -Mat4_Minor(*S, 2, 3)
			A\_24 =  Mat4_Minor(*S, 2, 4)
			A\_31 =  Mat4_Minor(*S, 3, 1)
			A\_32 = -Mat4_Minor(*S, 3, 2)
			A\_33 =   Mat4_Minor(*S, 3, 3)
			A\_34 = -Mat4_Minor(*S, 3, 4)
			A\_41 = -Mat4_Minor(*S, 4, 1)
			A\_42 =   Mat4_Minor(*S, 4, 2)
			A\_43 = -Mat4_Minor(*S, 4, 3)
			A\_44 =   Mat4_Minor(*S, 4, 4)
			det = Mat4_Det(A)
			k_det = 1/det
			Mat4_Scale(A, A, k_det)
			Mat4_Transpose(*T, A)
		EndProcedure
		Procedure		Mat4_CreateIdentity(*M.mat4)
			Mat4_SetL1(*M, 1, 0, 0, 0)
			Mat4_SetL2(*M, 0, 1, 0, 0)
			Mat4_SetL3(*M, 0, 0, 1, 0)
			Mat4_SetL4(*M, 0, 0, 0, 1)
		EndProcedure
		Procedure		Mat4_CreateRandom(*M.mat4)
			Mat4_SetL1(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
			Mat4_SetL2(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
			Mat4_SetL3(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
			Mat4_SetL4(*M, Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1), Real_Random(-1, 1))
		EndProcedure
		Procedure		Mat4_CreatePerspective(*M.mat4, fovy.f, aspect.f, zNear.f, zFar.f)
			Protected.f f, dz
			f = Cos(fovy/2)/Sin(fovy/2)
			dz = (zNear-zFar)
			Mat4_SetL1(*M, f/aspect, 0, 0, 0)
			Mat4_SetL2(*M, 0, f, 0, 0)
			Mat4_SetL3(*M, 0, 0, (zFar+zNear)/dz, 2*zFar*zNear/dz)
			Mat4_SetL4(*M, 0, 0, -1, 0)
		EndProcedure
		Procedure		Mat4_CreateOrtho(*M.mat4, left.f, right.f, bottom.f, top.f, zNear.f, zFar.f)
			Protected.f tx, ty, tz
			tx = -(right + left)/(right - left)
			ty = -(top + bottom)/(top - bottom)
			tz = -(zFar + zNear)/(zFar - zNear)
			Mat4_SetL1(*M, 2/(right - left), 0, 0, tx)
			Mat4_SetL2(*M, 0, 2/(top - bottom), 0, ty)
			Mat4_SetL3(*M, 0, 0, -2/(zFar - zNear), tz)
			Mat4_SetL4(*M, 0, 0, 0, 1)
		EndProcedure
		Procedure		Mat4_CreateTranslationComp(*M.mat4, vx.f, vy.f, vz.f)
			Mat4_SetL1(*M, 1, 0, 0, vx)
			Mat4_SetL2(*M, 0, 1, 0, vy)
			Mat4_SetL3(*M, 0, 0, 1, vz)
			Mat4_SetL4(*M, 0, 0, 0, 1)
		EndProcedure
		Procedure		Mat4_CreateTranslation(*M.mat4, *v.vec3)
			Mat4_CreateTranslationComp(*M, *v\x, *v\y, *v\z)
		EndProcedure
		Procedure		Mat4_CreateRotationComp(*M.mat3, rx.f, ry.f, rz.f)
			Protected.mat3 M3
			Mat3_CreateRotationComp(M3, rx, ry, rz)
			Mat4_SetM(*M, M3, Vec4_UnityW)
		EndProcedure
		Procedure		Mat4_CreateRotation(*M.mat3, *rot.vec3)
			Protected.mat3 M3
			Mat3_CreateRotation(M3, *rot)
			Mat4_SetM(*M, M3, Vec4_UnityW)
		EndProcedure
		Procedure		Mat4_CreateLookAt(*M.mat4, *source.vec3, *target.vec3, *up.vec3)
			Protected.vec3 dir, right, up2
			Protected.mat3 M3
			Protected.vec4 pos
			Vec3_Sub(dir, *target, *source)
			Vec3_Normalize(dir)
			Vec3_Cross(right, dir, *up)
			Vec3_Normalize(right)
			Vec3_Cross(up2, right, dir)
			Mat3_SetV(M3, right, up2, dir)
			Vec4_Vec(pos, *source, 1)
			Mat4_SetM(*M, M3, pos)
		EndProcedure
		Procedure		Mat4_FromQuat(*M.mat3, *Q.vec4)
			Protected.mat3 M3
			Mat3_FromQuat(M3, *Q)
			Mat4_SetM(*M, M3, Vec4_UnityW)
		EndProcedure
		Procedure		Mat4_Debug(*A.mat4)
			Debug StrF(*A\_11, 1) + ", " + StrF(*A\_12, 1) + ", " + StrF(*A\_13, 1) + ", " + StrF(*A\_14, 1)
			Debug StrF(*A\_21, 1) + ", " + StrF(*A\_22, 1) + ", " + StrF(*A\_23, 1) + ", " + StrF(*A\_24, 1)
			Debug StrF(*A\_31, 1) + ", " + StrF(*A\_32, 1) + ", " + StrF(*A\_33, 1) + ", " + StrF(*A\_34, 1)
			Debug StrF(*A\_41, 1) + ", " + StrF(*A\_42, 1) + ", " + StrF(*A\_43, 1) + ", " + StrF(*A\_44, 1)
		EndProcedure
	;}
;}
;{ Quaternions
	;{ variables
		Threaded _quat_mulv_mat.mat3
		Threaded _quat_divv_mat.mat3
		Threaded _quat_lerp_q1.vec4
		Threaded _quat_lerp_q2.vec4
		Threaded _quat_slerp_q1.vec4
		Threaded _quat_slerp_q2.vec4
		Threaded _quat_slerp_w.f
		Threaded _quat_slerp_inv_sin_w.f
	;}
	;{ Procedures
		Macro			Quat_Conjugate(Qr, Qa)
		    Qr\x =  -Qa\x
		    Qr\y = -Qa\y
		    Qr\z = -Qa\z
		    Qr\w = Qa\w
		EndMacro
		Procedure	Quat_FromAxisComp(*Q.vec4, ax.f, ay.f, az.f, angle.f)
			Protected.f sin_a, cos_a
			sin_a = Sin(angle/2)
			cos_a = Cos(angle/2)
			*Q\x = ax * sin_a
			*Q\y = ay * sin_a
			*Q\z = az * sin_a
			*Q\w = cos_a
		EndProcedure
		Procedure	Quat_FromAxis(*Q.vec4, *axis.vec3, angle.f)
			Quat_FromAxisComp(*Q, *axis\x, *axis\y, *axis\z, angle)
		EndProcedure
		Procedure	Quat_FromMat(*Q.vec4, *M.mat3)
			Protected.f tr
			tr = *M\_11 + *M\_22 + *M\_33
			If tr > 0
				Vec4_Set(*Q, *M\_23 - *M\_32, *M\_31 - *M\_13, *M\_12 - *M\_21, tr + 1)
			Else
				If *M\_11 > *M\_22 And *M\_11 > *M\_33
					Vec4_Set(*Q, 1 + *M\_11 - *M\_22 - *M\_33, *M\_21 + *M\_12, *M\_31 + *M\_13, *M\_23 - *M\_32)
				Else
					If *M\_22 > *M\_33
						Vec4_Set(*Q, *M\_21 + *M\_12, 1 + *M\_22 - *M\_11 - *M\_33, *M\_32 + *M\_23, *M\_31 - *M\_13)
					Else
						Vec4_Set(*Q, *M\_31 + *M\_13, *M\_32 + *M\_23, 1 + *M\_33 - *M\_11 - *M\_22, *M\_12 - *M\_21)
					EndIf
				EndIf
			EndIf
			Vec4_Normalize(*Q)
		EndProcedure
		Procedure	Quat_MulQ(*Qr.vec4, *Q1.vec4, *Q2.vec4)
 			Vec4_Hamilton(*Qr, *Q1, *Q2)
 		EndProcedure
		Macro			Quat_MulV_Setup(Q)
			Mat3_FromQuat(_quat_mulv_mat, Q)
		EndMacro
		Macro			Quat_MulV(Vr, V)
			Mat3_MulV(Vr, _quat_mulv_mat, V)
		EndMacro
		Macro			Quat_DivV_Setup(Q)
			Mat3_FromQuat(_quat_mulv_mat, Q)
			Mat3_Transpose(_quat_divv_mat, _quat_mulv_mat)
		EndMacro
		Macro 			Quat_DivV(Vr, V)
			Mat3_MulV(Vr, _quat_divv_mat, V)
		EndMacro
		Procedure	Quat_Lerp_Setup(*Q1.vec4, *Q2.vec4)
			Protected.f inner
			Vec4_Copy(_quat_lerp_q1, *Q1)
			Vec4_Copy(_quat_lerp_q2, *Q2)
			inner = Vec4_Dot(_quat_lerp_q1, _quat_lerp_q2)
			If inner < 0
				Vec4_Scale(_quat_lerp_q2, _quat_lerp_q2, -1)
			EndIf
			Vec4_Sub(_quat_lerp_q2, _quat_lerp_q2, _quat_lerp_q1)
		EndProcedure
		Macro			Quat_Lerp(Q, t)
	    	Vec4_Copy(Q, _quat_lerp_q2)
	    	Vec4_Scale(Q, t)
	    	Vec4_Add(Q, Q, _quat_lerp_q1)
		EndMacro
		Procedure	Quat_Slerp_Setup(*Q1.vec4, *Q2.vec4)
			Protected.f cos_w
			Vec4_Copy(_quat_slerp_q1, *Q1)
			Vec4_Copy(_quat_slerp_q2, *Q2)
			cos_w = Vec4_Dot(_quat_slerp_q1, _quat_slerp_q2)
			If cos_w < 0
				cos_w = -cos_w
				Vec4_Scale(_quat_slerp_q2, _quat_slerp_q2, -1)
			EndIf
			If cos_w > 0.9999
				cos_w = 0.9999
			EndIf
			_quat_slerp_w = ACos(cos_w)
			_quat_slerp_inv_sin_w = 1/Sin(_quat_slerp_w)
			Vec4_Scale(_quat_slerp_q1, _quat_slerp_q1, _quat_slerp_inv_sin_w)
			Vec4_Scale(_quat_slerp_q2, _quat_slerp_q2, _quat_slerp_inv_sin_w)
		EndProcedure
		Procedure	Quat_Slerp(*Q.vec4, t.f)
			Protected.vec4 q1, q2
			Protected.f k
			k = Sin((1-t)*_quat_slerp_w)
			Vec4_Scale(q1, _quat_slerp_q1, k)
			k = Sin(t*_quat_slerp_w)
			Vec4_Scale(q2, _quat_slerp_q2, k)
			Vec4_Add(*Q, q1, q2)
		EndProcedure
	;}
;}



DisableExplicit






























; IDE Options = PureBasic 5.11 (Windows - x86)
; CursorPosition = 178
; Folding = CAIAABAEAACEA+AAAAAAAAAAgAAg+
; EnableXP